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We report results of quantum Monte Carlo simulations in the canonical and the grand-canonical 
ensemble of the two- and three-dimensional Bose-Hubbard model with quadratic and quartic con- 
fining potentials. The quantum criticality of the superfluid-Mott insulator transition is investigated 
both on the boundary layer separating the two coexisting phases and at the center of the traps where 
the Mott-insulating phase is first established. Recent simulations of systems in quadratic traps have 
shown that the transition is not in the critical regime due to the finite gradient of the confining 
potential and that critical fluctuations are suppressed. In addition, it has been shown that quantum 
critical behavior is recovered in flat confining potentials as they approach the uniform regime. Our 
results show that quartic traps display a behavior similar to quadratic ones, yet locally at the center 
of the traps the bulk transition has enhanced critical fluctuations in comparison to the quadratic 
case. Therefore quartic traps provide a better prerequisite for the experimental observation of true 
quantum criticality of ultracold bosonic atoms in optical lattices. 

PACS numbers: 03.75.Lm, 73.43.Nq, 05.30.Jp 



I. INTRODUCTION 



The Bose-Hubbard model offers an almost perfect de- 
scription of Bose-Einstein condensates in optical lat- 
tices 0, 0. Therefore the model has been intensively 
studied both analytically El H, El H| as well as numer- 
ically [1 H H [H EUTIl El hi recent years. In 
the absence of a trapping potential, i.e., for the homoge- 
neous Bose-Hubbard model, a quantum phase transition 
from a superffuid to a Mott-insulating phase occurs at 
commensurate fillings upon increasing the optical lattice 
depth Q . Experimentally, Bose-Einstein condensates are 
created by cooling bosonic atoms. In order to prevent 
them from evaporating, they must be trapped by a con- 
fining potential. As shown for example in Ref. fill inho- 
mogeneous confining potentials induce a coexistence of 
superfluid and Mott-insulating regions for a continuous 
range of incommensurate fillings. By a suitable increase 
of the lattice depth and the number of particles, one may 
start forming a Mott-insulating region in the center of 
the trap. A further increase of the lattice depth induces 
a change of the volume fractions between the two phases. 
Accordingly, the transition recently observed in experi- 
ments of ultracold Bose gases in optical lattices embed- 
ded in quadratic confining potentials 0] should be better 
viewed as a crossover with changing volume fractions of 
the two phases [l2j rather than as a true phase transition. 

Although globally no true transition can be observed 
in systems with inhomogeneous traps, one might expect 
to observe a superfluid-Mott insulator transition on the 
boundary layer separating the coexisting superfluid and 



Mott-insulating regions and, locally, at the center of the 
trap where the Mott-insulating phase has been estab- 
lished. In the following we distinguish these transitions 
by referring to them as "surface transition" and "bulk 
transition," respectively. Recent quantum Monte Carlo 
simulations of the Bose-Hubbard model in quadratic con- 
fining potentials have shown that the quantum criticality 
of these transitions is destroyed by the finite gradient of 
the confining potential ^| • However it has been demon- 
strated that quantum criticality is recovered in flat traps 
where the potential gradient becomes irrelevant in the 
center, i.e., when locally the trap center approaches the 
uniform regime. 

Here we study flat anharmonic traps, which can be 
experimentally realized, e.g., by superimposing pairs of 
weak, repulsive, blue-detuned Gaussian laser beams to 
pairs of attractive, red-detuned Gaussian laser beams 
that create the optical lattice. Depending on the beam 
parameters, it is possible to cancel out the harmonic 
terms in the series expansion of the resulting potential. 
In our simulations with anharmonic traps we consider 
only fourth-order terms in the expansion and neglect all 
higher-order ones. The results we have obtained from 
our simulations of two-dimensional (2D) systems in quar- 
tic traps and the comparison to our results for systems 
in quadratic traps, closed-box systems confined by com- 
pletely flat and infinitely sharp traps, as well as ho- 
mogeneous systems with periodic boundary conditions 
clearly demonstrate that the surface transition in sys- 
tems with quartic traps has the same lack of quantum 
criticality as the corresponding transition in the sys- 
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terns with quadratic traps. However, the bulk transition 
has stronger critical fluctuations in the center of quartic 
traps, than in the center of quadratic traps. These obser- 
vations lead us to the conclusion that quantum critical- 
ity in systems with inhomogencous traps is determined 
by the flatness in the trap centers and their closest sur- 
roundings, i.e., by the smaller deviation from the uniform 
regime, and not by the shape and steepness of the walls. 
Thus, quartic traps could provide a better prerequisite for 
future experimental observations of true quantum criti- 
cality. 

We also present results on three-dimensional (3D) sys- 
tems with smaller system sizes which qualitatively show 
the same results as in two space dimensions, albeit less 
pronounced. 

The paper is organized as follows: In Sec. ^ we in- 
troduce in detail both the Bose-Hubbard model with 
the specific traps that we have used and the observables 
needed for the investigations of the systems. In Sec. IIIII 
we present our results, and in Sec. II VI we summarize our 
observations. 



II. MODEL, OBSERVABLES, AND 
NUMERICAL DETAILS 



In order to describe trapped ultracold bosonic atoms 
in optical lattices we use the Bose-Hubbard Hamiltonian 

n = 1 E ( s & + H - c -) + 1 E - x ) 



(id) 

MoE^* + E V(i)ni> 



(1) 



where b\ and hi are the creation and annihilation opera- 
tors for bosons at lattice sites i, respectively, and 



bib, 



(2) 



is the local density operator, t represents the nearest- 
neighbor hopping matrix element, U the on-site repul- 
sion, no a chemical potential offset that controls the fill- 
ing of the trap and V(i) the trap potential. The two last 
terms of the Hamiltonian can be conveniently combined 
by defining an effective, spatially-dependent chemical po- 
tential fif that is experienced by a boson at site i, i.e., 



;ff = w - v(i) 



(3) 



Due to the inhomogeneity introduced by a trapping 
potential needed for the confinement of the atoms in 
experimental realizations of the Bose-Hubbard model, 
it becomes necessary to consider local, site-dependent 
quantities, in contrast to the homogeneous Bose-Hubbard 
model, where it suffices to measure global quantities 
such as the global compressibility and the superfluid den- 
sity 01 to capture the system characteristics. One way 
of distinguishing coexisting local superfluid and Mott- 
insulating states in systems with traps is to investigate 
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FIG. 1: (Color online) Evolution of the measured particle 
number N (solid curve) toward the desired value N = 850 
while adjusting the chemical potential offset fio/U (dashed 
curve) at the beginning of the thermalization phase. The 
system size is 50 x 50 in a quartic trap. The inset is a detailed 
view of N fluctuating around N = 850 as a function of time. 
The circles indicate the configurations for which N is exactly 
850 and which can be used in a canonical measurement. 



the topology of local density profiles. Regions with con- 
stant integer fillings can be interpreted to be Mott insu- 
lating, whereas regions with noninteger fillings must be 
superfluid. The local compressibility K' ocal at site i 



„ local 



d(N) 



elf 



&T[{ni{T)N)-{ni{T)){N)], (4) 



provides a more precise way to distinguish the local states 
and to probe for quantum criticality. It expresses the 
response of the system's particle number N to a local 
change of the effective chemical potential fi^ s at site i. 



In the above equation, /3 
temperature 



1/ksT denotes the inverse 



hi{r) = exp(r7i)nj exp(— tTL)\ 



(5) 



the imaginary-time propagated operator, and (• • •} the 
Monte Carlo sample average. Note that this definition of 
the local compressibility only makes sense in the grand- 
canonical ensemble where the total number of particles 
N in the system is variable, since otherwise the integrand 
vanishes and the local compressibility is zero everywhere. 

The measurement of local density fluctuations provides 
another way of investigating the local states of inhomo- 
geneous systems. In addition to measuring the variance 
of the local density 



A, 



K 2 ) - <^> 2 



(6) 



it can be useful to measure the on-site compressibility 
K° nslte , i.e., the response of the local density at site i to 
a chemical potential change at this site 



djrii) 
d/jf 







dT[(f»i(r)ni(0)> - (n i (r))(n,(0))] 



(7) 

For our simulations, we have used two wormlike quan- 
tum Monte Carlo algorithms provided by the ALPS 
project |15| , the worm algorithm jl6j | , and the stochastic 
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TABLE I: Overview of the simulations with quadratic and 
quartic traps in 2D simulated in the canonical ensemble. N is 
the fixed particle number at which the observables are evalu- 
ated in the canonical ensemble, (j,q/U and t/U are the chem- 
ical potential offset and the hopping parameter in units of 
the interaction parameter U, respectively. a 2 and 04 are the 
curvatures of the quadratic and quartic traps. 



Size 


Trap curvature 


N 


fto/U 


t/U 


50 2 


ai/U = 3.5 ■ 1CT 6 


1200 


0.362 


0.04 


50 2 


a 2 /U = 2 ■ HT 3 


600 


0.370 


0.04 



series expansion (SSE) algorithm with directed loops [T^| . 
Both algorithms work gencrically in the grand-canonical 
ensemble. Since experiments are usually performed with 
a fixed number of particles, we want our simulations to 
also be performed with a fixed particle number, i.e., in the 
canonical ensemble. For this purpose, we have applied 
some modifications to the algorithms to fit our needs. 
In order for the Hamiltonian in Eq. to be fully de- 
fined, the correct chemical potential offset fio/U, which 
yields the desired particle number TV at which the simula- 
tions are to be performed, has to be determined. Starting 
from a guess value the correct value of /j.q/U is approx- 
imated at the beginning of the thermalization phase of 
the simulation, as illustrated in Fig.^ The dashed curve 
shows the subsequent adjustments of the chemical poten- 
tial offset no/U and the solid curve represents the evolu- 
tion of the measured particle number, which equilibrates 
around the desired particle number N (solid, horizontal 
line) once the adjustment of /iq/U is complete. 

Measuring observables in the canonical ensemble is 
then simply done by measuring the particle number after 
each Monte Carlo step and by selecting out only updated 
configurations for which the measured particle number 
coincides with the desired particle number at which the 
simulations are meant to be performed. By computing 
the Monte Carlo sample averages using only the selected 
updated configurations, the observables are evaluated in 
the canonical ensemble, whereas by including all other 
updated configurations for the computation of the sam- 
ple averages, the observables are evaluated in the grand- 
canonical ensemble. 

We have found no substantial qualitative deviations 
between the canonical and grand-canonical data for all 
observables. In Figs. [3| and 0] we plot the data evaluated 
in the canonical ensemble, except for the local compress- 
ibility, which, for reasons mentioned above, cannot be 
measured canonically. In Figs. [S] [7| and |S| where com- 
parisons are made with uniform systems, we have chosen 
to plot the data evaluated in the grand-canonical ensem- 
ble, since obtaining the different curves for the uniform 
systems implies that the filling is allowed to change . 

The 2D and 3D trapping potentials we have used are 



TABLE II: Overview of the simulations with quadratic and 
quartic traps in 2D (top) and 3D (bottom) simulated in the 
grand-canonical ensemble. TV is the particle number and 
Mcfnter/t^ an d t/U are the effective chemical potential at the 
central lattice sites [l9|| and the hopping parameter in units 
of the interaction parameter U, respectively. 0,2 and 04 are 
the curvatures of the quadratic and quartic traps. 



Size 


Trap curvature 


N 


/^center i^J 


t/U 


50 2 
50 2 


cla/U = 3.5 ■ 10" 6 
a 2 /U = 2- 10" 3 


908 - 924 
341 - 413 


0.217 
0.217 


0.01 - 0.06 
0.01 - 0.06 


14 3 
14 3 


m/U = 7 ■ 10 -4 
a 2 /U = 3- 10" 2 


462 - 498 
90 - 122 


0.217 
0.217 


0.006 - 0.04 
0.006 - 0.04 


1.5 
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FIG. 2: (Color online) Sections V(x,y = a/2) and V(x,y = 
a/2, z = a/2) through the 2D (left) and 3D (right) traps that 
have been used in the simulations. The simulation parameters 
are summarized in Tables and EI The points indicate the 
lattice sites separated by the lattice constant a. 

V$ = a 2 (x 2 + y 2 + z 2 ) , (10) 
V$ = a 4 (x 4 + y 4 + z 4 ), (11) 

where the trap curvatures are given in Tables [I] and [n] 
(see also Fig. [2J. The lattice sites in the 2D and 3D 
systems in units of the lattice spacing a are located at 
the coordinates (x,y) — (±[n + l/2],±[n+ 1/2]) and 
(x,y,z) = (±[n+ 1/2], ±[71 + l/2],±[n + 1/2]), respec- 
tively, where n — 0, l,...,L/2— 1. The linear ex- 
tents of the systems with quadratic and quartic traps 
are L/a = 50 in 2D and L/a = 14 in 3D. With the trap 
curvatures and the range of the chemical potential offsets 
tM)/U that we have used, these system sizes ensure that 
the effective chemical potential takes large enough nega- 
tive values at the boundaries of the systems for the par- 
ticles to never reach the boundaries. The linear extents 
of the systems with closed boxes and the homogeneous 
systems with periodic boundary conditions have been re- 
duced to L/a = 32 and L/a = 28 in 2D and L/a = 8 
in 3D to approximatively match the effective number of 
bosons trapped in the quadratic and quartic traps, i.e., 
the sizes of the regions with a nonvanishing local density. 



III. RESULTS 



= a 2 (x 2 + y 2 ) , 
= a 4 (x 4 + y A ) , 



(8) 
(9) 



In this section we present the results obtained from 
simulations of the two- and three-dimensional Bose- 
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FIG. 3: (Color online) Spatial dependence of the local den- 
sity n measured in the canonical ensemble (top), and of the 
local compressibility K local f measured in the grand-canonical 
ensemble (bottom) for 2D systems of size 50 x 50 containing 
respectively 600 bosons in a quadratic trap with curvature 
a-z/U — 2 ■ 10 -3 , chemical potential offset no/U = 0.370 and 
hopping parameter t/U = 0.04 (left) and 1200 bosons in a 
quartic trap with curvature cn/U = 3.5 • 10~ 6 , chemical po- 
tential offset no/U = 0.362 and the same hopping parameter 
t/U = 0.04 (right). The special fourfold symmetry of the 
quartic trap is clearly reflected by the shape of the superfluid 
shell between the outer n — and the central n — 1 Mott 
plateaus (right). 



Hubbard model in quartic traps and compare results for 
quadratic traps with closed box configurations and ho- 
mogeneous systems with periodic boundary conditions. 
Tables [|] and [H] give an overview of the simulation pa- 
rameters in the canonical and grand-canonical ensemble, 
respectively. 

First, we investigate quantum criticality at the sur- 
face which separates the superfluid and Mott-insulating 
states, which is shell structured in systems with quadratic 
and quartic traps. Our results suggest that in regards to 
quantum criticality at the surface quadratic and quartic 
traps are equivalent: both systems show the same lack 
of quantum criticality. Next, we study quantum critical- 
ity to the center of the traps by restricting the measure- 
ments to small regions around the center. We drive these 
regions through the superfluid-Mott insulator transition 
and show that this bulk transition in the center of quartic 
traps has enhanced critical fluctuations in comparison to 
the center of quadratic traps. 



A. Quantum criticality on the boundary between 
superfluid and Mott-insulating regions 

We compare a system with a quadratic trap contain- 
ing 600 bosons and a system with a quartic trap contain- 
ing 1200 bosons. The hopping parameter is t/U = 0.04 
in both cases and the chemical potential offsets are 
Ho/U — 0.370 in the quadratic case and fio/U — 0.362 in 



the quartic case. In both cases a Mott-insulating plateau 
with integer density (rii) = 1 is present at the center of 
the trap, which is surrounded by a superfluid ringlike re- 
gion, as depicted in Fig. [21 where the spatial dependence 
of the local density n and of the local compressibility 
k , is shown. The left panels show the quadratic case 
and the right panels the quartic case. In both cases, the 
local compressibility takes its largest values close to the 
outer boundary of the superfluid shell. The shape of the 
superfluid shell clearly reflects the special fourfold rota- 
tional symmetry of the quartic trap, compared to the 
quadratic trap, which has a continuous rotational sym- 
metry. Furthermore, the width of the superfluid shell is 
clearly smaller in the quartic trap than in the quadratic 
trap, which is due to the larger steepness of the former 
in the corresponding regions of the superfluid shell. 

A better quantitative description of the systems in the 
two different traps is given by the radial dependence of 
the observables from the center towards the boundaries 
of the systems. In the quadratic case, due to the contin- 
uous rotational symmetry, the values at all lattice sites 
can be used to create such geometric profiles. In the 
quartic case, we are limited by the fourfold symmetry 
of the underlying potential and we thus use only points 
lying on the sections through the center and parallel to 
the boundaries of the systems. This limitation no longer 
applies to quartic traps when we plot local quantities as a 
function of the local effective chemical potential /i eS /U, 
since, as follows from the data collapse on single curves, 
cf. Fig. [S] a local potential approximation holds just as 
for quadratic traps [lj], i.e., these quantities can be de- 
termined from the value of the local effective chemical 
potential. 

The errors have been determined by means of a non- 
parametric bootstrap analysis [2(j applied to each set of 
values corresponding to the same radial distance from 
the center or to the same effective chemical potential. If 
the error bars are not visible, they are smaller than the 
corresponding point sizes. 

The radial dependence of the local density n, the local 
compressibility K local , the on-site compressibility K onslte 
and the variance A for the same systems introduced 
above are given in Fig.0] Noninteger local densities cor- 
respond to regions where the system is in the superfluid 
state, whereas integer local densities indicate regions in 
which the system is in the Mott-insulating state. The 
local compressibility profile consists of an asymmetric 
double peak, which reflects the increase of the particle 
number fluctuations near the boundaries of the Mott- 
insulating regions. Both peaks are of the same height in 
the corresponding hard-core model due to particle-hole 
symmetry |l2]. The on-site compressibility, K onsltc , and 
the variance, A, both peak inside the superfluid shell, but 
in contrast to the local compressibility, K local , they do not 
completely vanish inside the central Mott plateau. This 
is due to virtual hopping processes that are completely 
suppressed only by an infinite energy gap in the limit 
t/[/-»O0. As for the quadratic traps [3, 

we find 
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FIG. 4: (Color online) Radial dependence of the local density 
n, the local compressibility K local , the on-site compressibility 
^.onsitc anc j ^ ne var j ance /\ f or ^ ne same systems as in Fig. 

All quantities are measured in the canonical ensemble, except 
for K local that can only be measured in the grand-canonical 
ensemble. 



for the quartic traps that K local is a better probe than 
^onsite an( j ^ £ Qr existence of superfiuid and Mott- 
insulating regions. It can, therefore, serve as a genuine 
order parameter to characterize the superfiuid-Mott in- 
sulator transition. 

In Fig.[S]we show the local density n and the local com- 
pressibility functions of the local effective chemi- 
cal potential /z eff /U, which, as pointed out in Ref.0, are 
not universal functions, but depend on both the geome- 
tries and the curvatures of the traps. The data for both 
quantities taken from the simulations with quadratic and 
quartic traps match almost perfectly. As a comparison, 
both curves are plotted for closed-box systems confined 
by completely flat and infinitely sharp traps as well as 
for homogeneous systems with periodic boundary condi- 
tions, all of size 32 x 32. This smaller size is comparable to 
the effective sizes of the systems with the quadratic and 
the quartic traps. Since in these uniform cases the effec- 
tive chemical potential remains constant over the whole 
system, a whole set of simulations, each of which is per- 
formed with a different chemical potential, is needed in 
order to obtain these additional curves. In the local den- 
sity curve, cusps appear at the points where the den- 
sity approaches n — and n — 1 and correspondingly 
both singularities in the local compressibility are more 
pronounced, indicating quantum criticality. As shown in 
RefM 

these cusps in n(fj, cS /U) become smoother as the 
curvature is increased from zero (uniform case) to finite, 
positive values. This is due to the gradient in the con- 
fining potential that becomes increasingly relevant and 
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FIG. 5: (Color online) Local density n (top) and local com- 
pressibility K local (bottom), measured in the grand-canonical 
ensemble, vs the local effective chemical potential /i off /V . The 
quadratic and quartic trap cases are the same as Fig. [3] We 
also show systems in closed boxes of size 32 x 32 as well 
as homogeneous systems with periodic boundary conditions 
also of size 32 x 32. In all system, the hopping parameter is 
t/U = 0.04. 



which destroys the quantum criticality. No quantum crit- 
ical behavior is found for the transition at the surface 
separating the coexisting superfiuid and Mott-insulating 
states in the quartic trap, in spite of its flatness in the 
center and the larger steepness of its walls. 



B. Bulk quantum criticality in the trap center 

Instead of investigating the transition at the surface 
separating the coexisting superfiuid and Mott-insulating 
states in the systems with quadratic and quartic traps, we 
now focus on the center of these traps and drive the local 
state through the superfiuid-Mott insulator transition, 
which can be induced by either increasing the chemical 
potential offset (j,q/U or by decreasing the ratio t/U be- 
tween the hopping and the interaction parameter. The 
latter case, which is more likely to be reproduced experi- 
mentally, is shown in Fig.0 where the local density n and 
the local compressibility K local are plotted against t/U at 
constant chemical potential fi/U — 0.217. This scan is 
marked with a dashed line in Fig. where the phase di- 
agram of the homogeneous Bose-Hubbard model Q is 
depicted schematically. At each point of the scan, both 
n and K local are measured at the twelve lattice sites that 
are closest to the center, i.e., at the four central lattice 
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FIG. 6: (Color online) Phase diagram of the homogeneous 
Bose-Hubbard model |3j. It consists of a series of Mott- 
insulating lobes that are surrounded by a superfluid region. 
Each of the Mott lobes is characterized by an integer filling of 
single lattice sites. Here, only the n = and n — 1 Mott lobes 
are shown. The constant (fl/U — 0.217) scans presented in 
Figs. and [H] are indicated by the dashed, horizontal line. 



sites pjj and at their eight nearest neighbors (cross con- 
figuration). In 2D, with our chosen system sizes and 
trap curvatures, the difference of the effective chemical 
potential felt by bosons on the four central lattice sites 
and the one felt by bosons on their nearest neighbor lat- 
tice sites is negligible and thus the values measured at 
the nearest neighbor lattice sites can be included to en- 
hance the statistics without introducing any bias. The 
curves originating from the system with the quartic trap 
almost match the curves originating from the uniform 
systems with closed boxes and periodic boundary con- 
ditions, whereas for the curves obtained from the sys- 
tem with the quadratic trap, the deviations are much 
stronger. The sharp cusp in the K local curve of the sys- 
tem with the quartic trap indicates the increase of critical 
fluctuations inside the trap and will diverge in the con- 
tinuum limit. This special behavior demonstrates that 
locally at the center, the bulk transition from the super- 
fluid to the Mott-insulating state exhibits stronger criti- 
cal fluctuations in the quartic trap than in the quadratic 
trap. 

Similarly as in 2D, we have performed the same con- 
stant (n/U — 0.217) scans through the critical point 
also for 3D systems, see Fig. [(J In contrast to the 2D 
case, the local density n and the local compressibility 
«.iocai have been measured only at the eight central lat- 
tice sites Due to the smaller system sizes and the 
larger trap curvatures than in 2D, the twenty- four near- 
est neighbors lattices sites could not be considered since 
the effective chemical potentials at these nearest neigh- 
bors lattice sites and at the central lattice sites differ 
significantly from each other and thus correspond to two 
different scan lines in the phase diagram. As can be ex- 
tracted from Fig. |S| the same yet less pronounced be- 
havior in regards of bulk quantum criticality is observed, 
i.e., locally at the center of the traps, the transition from 
the local superfluid to the local Mott-insulating state is 
found to have slightly stronger critical fluctuations in the 
quartic than in the quadratic trap. However, care must 
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FIG. 7: (Color online) Constant (fi/U = 0.217) scans through 
the critical point. The density n (top) and the local com- 
pressibility K local (bottom), in the grand-canonical ensemble, 
vs t/U for 2D systems of 50 x 50 sites. The quadratic trap 
curvature is (12/U = 2 • 10 -3 and the quartic trap curvature 
ai/U = 3.5 ■ 10 -6 . Both quantities are measured at the cen- 
ter of the traps. For comparison, the same scans are plotted 
for a system in a closed-box and one with periodic boundary 
conditions of size 28 x 28. 



be taken while interpreting these results, due to the fi- 
nite size effects that strongly affect the centers of these 
small 3D systems with their significantly reduced linear 
extents, compared to the 2D systems presented above. 



IV. CONCLUSIONS 

Ultracold bosonic atoms in optical lattices have been 
studied in quadratic and quartic traps. The physics of 
these confined systems is almost perfectly captured by 
the Bose-Hubbard model [jj. The specific traps can be 
taken into account by a site-dependent effective chemical 
potential in the Bose-Hubbard Hamiltonian. Our results 
in 2D clearly demonstrate that the quantum criticality 
of the transition at the surface layer separating the coex- 
isting superfluid and Mott-insulating states is destroyed 
by the finite gradient of the trapping potential, even for 
traps with quartic profiles, which have steep walls and in 
whose centers the gradient almost vanishes. However, we 
have found that, if the measurements are restricted to the 
center of the traps and the local state is driven through 
the superfluid-Mott insulator transition, the then ob- 
served bulk transition has stronger critical fluctuations 
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FIG. 8: (Color online) Constant (n/U = 0.217) scans through 
the critical point. The density n (top) and the local com- 
pressibility K local (bottom), in the grand-canonical ensemble, 
are plotted for 3D systems of 14 3 sites in a quadratic trap 
(curvature (12/U = 0.03) and in a quartic trap (curvature 
a 4 /U = 7 ■ 10~ 4 ). As in the 2D scans (Fig.0, n and K local 
are measured at the center of the traps. For comparison, the 
same scans are plotted for closed-box systems and homoge- 
neous systems with periodic boundary conditions of size 8 3 , 
which approximatively corresponds to the effective sizes of the 
systems with the quadratic and quartic traps. 



in quartic traps than in quadratic traps. To demonstrate 
this, we use the local compressibility as a genuine, lo- 



cal order parameter to characterize the transition and 
drive the system center from the superfluid to the Mott- 
insulating state by decreasing the ratio t/U between the 
hopping and the interaction parameter at constant chem- 
ical potential /i/U '. The deviations of the curves obtained 
from the systems with quartic traps from the curves ob- 
tained from uniform systems with closed boxes and peri- 
odic boundary conditions are much less pronounced than 
the deviations of the curves obtained from the systems 
with quadratic traps. Even if there is no perfect match, 
we expect the quartic traps to constitute a remarkably 
better prerequisite for the experimental detection of true 
quantum criticality. 

Our simulations of moderately-sized 3D systems lead 
to the same conclusions as in 2D, but, due to strong finite 
size effects present in these systems, an analysis with yet 
larger systems is desirable. 

Finally, it is worth mentioning that quartic traps rep- 
resent an ideal for traps that are experimentally realized 
by superimposing blue- to red-detuned laser beams with 
Gaussian intensity profiles. The flatness achievable in 
the trap center depends on the ability to suppress the 
quadratic terms in the series expansion of the resulting 
confining potential. Since this might turn out to be a 
rather difficult task, it could be interesting to investigate 
what influence small fluctuations in the otherwise flat 
trap centers have on the local bulk quantum criticality. 
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